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We study hysteresis for a two-dimensional, spin-1/2, nearest-neighbor, kinetic Ising ferromagnet 
in an oscillating field, using Monte Carlo simulations and analytical theory. Attention is focused on 
^ , small systems and weak field amplitudes at a temperature below Tc- For these restricted parameters, 

■ the magnetization switches through random nucleation of a single droplet of spins aligned with the 

applied field. We analyze the stochastic hysteresis observed in this parameter regime, using time- 
dependent nucleation theory and the theory of variable-rate Markov processes. The theory enables 
us to accurately predict the results of extensive Monte Carlo simulations, without the use of any 
adjustable parameters. The stochastic response is qualitatively different from what is observed, 
either in mean-field models or in simulations of larger spatially extended systems. We consider the 
^ ' frequency dependence of the probability density for the hysteresis-loop area and show that its average 

slowly crosses over to a logarithmic decay with frequency and amplitude for asymptotically low 
frequencies. Both the average loop area and the residence-time distributions for the magnetization 
show evidence of stochastic resonance. We also demonstrate a connection between the residence- 
time distributions and the power spectral densities of the magnetization time series. In addition to 
their significance for the interpretation of recent experiments in condensed-matter physics, including 
studies of switching in ferromagnetic and ferroelectric nanoparticles and ultrathin films, our results 
are relevant to the general theory of periodically driven arrays of coupled, bistable systems with 
stochastic noise. 
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CN . I. INTRODUCTION 

> 

' Hysteretic response to an oscillating control parameter or "force" is a nonlinear nonequilibrium phenomenon com- 
I monly observed in both natural and man-made systems. The example most familiar in physics and electrical engi- 
fS| ' neering is probably the hysteresis loop produced when a ferromagnet at a temperature below its critical temperature 
^-H Tc is placed in an oscillating magnetic field . Similar behavior is seen in ferroelectrics . Some other examples 
are electrochemical adsorbate layers that are driven through a phase transition by an oscillating electrode potential in 
a Cyclic Voltammetry experime nt P ,|9| , systems driven through a phase transition between different liquid-crystalline 
phases by pressure oscillations ]10(|, and systems driven through a solid-liquid phase transition by temperature os- 
cillations. Hysteresis is often modeled by systems of differential equations that display discontinuous bifurcations 

Systems that exhibit hysteresis have in common a nonlinear, irreversible response, which causes the phase of the 
(~| response to lag behind the force. The physical mechanism that causes the hysteretic behavior can, however, be 
O 1 quite different in different systems and even in different parameter regimes for the same system. The details of this 
J-^ ■ mechanism must be considered in order to accurately predict such aspects of the hysteretic response as its dependence 
^ , on the frequency and amplitude of the oscillating force. Here we present a study of hysteresis in a particular model 
k>( ' system which incorporates both spatial degrees of freedom and thermal fluctuations - a kinetic Ising ferromagnet - 
rS in a parameter regime where the model has a first-order phase transition in equilibrium and the system response is 
stochastic. The model and its behavior in this regime are relevant to at least two different research areas that are 
■ ■ ' rarely discussed together: experimental studies of switching dynamics in nanoscale ferromagnetic and ferroelectric 
particles and ultrathin films, and theoretical and experimental studies of stochastic resonance in spatially extended 
systems. We hope the present study may contribute to some intellectual cross-fertilization. 

In recent years new experimental techniques, such as magnetic force microscopy (MFM) [^-^, have been de- 
veloped that permit measurements of the magnetization state and switching behavior of particles as small as a few 
nanometers. Ferromagnetic particles in this size range consist of a single domain in equilibrium, and together with 
ultrathin films they are of interest as potential materials for ultra-high density recording media. The dynamics of 
magnetization reversal in such systems has been modeled with kinetic Ising systems subject to sudden field reversal 
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|[19|-p3[ . These numerical and analytical studies have given results in qualitative agreement with the experiments men- 
tioned above. Recent experiments on ultrathin ferromagnetic Fe/Au (001) films [|4| have considered the frequency 
dependence of hysteresis loop areas, which were interpreted in terms of effective exponents consistent with those found 
for a continuous spin model [|7| j25| - [2^ . Similar experiments on ultrathin Co films on Cu(OOl) have found exponents 
consistent with a mean- field treatment of the Ising model ||2^. These studies of nanoparticles and ultrathin films 
suggest that experiments can now be performed on systems sufficiently small that atomic-scale simulations become 
feasible, and that kinetic Ising systems are useful models for switching in such nanoscopic systems. 

Since its introduction as a possible model for the time dependence of the Earth's ice ages [Q, the concept of 
stochastic resonance has been applied to a variety of phenomena in physical and biological science and engineering. 



in which response to a periodic force is enhanced by noise |30|. Most early treatments considered a single bistable 
element similar to a mean- field model of a ferromagnet plHSS" 



with added noise. However, more recently experimental 
studies have been conducted with chains of coupled diode resonators , and numerical and theoretical studies have 
considered locally coupled one-dimensional time-dependent Ginzburg-Landau or Frenkel-Kontorova models [p5|-p8|, 
Ising models in one two |^^, and three dimensions, chains of coupled nonlinear maps [Q, and systems of 
globally coupled bistable elements 

Here we consider hysteresis in a two-dimensional spin-1/2, nearest-neighbor, kinetic Ising ferromagnet in an oscil- 
lating field with periodic boundary conditions. For convenience, and because of the many experimental measurements 
of hysteresis that address magnetic systems, we use the customary magnetic language, in which the order parameter 
is the magnetization per site, m{t) S [—1,-1-1], and the force is the magnetic field H{t). However, we expect our 
results also to apply to stochastic hysteresis phenomena in other areas of science. For example, in dielectrics m(t) and 
H{t) can be re-interpreted as polarization and electric field, in adsorption problems as coverage d{t) = [2m{t) — 1] 
and (electro) chemical potential or (osmotic)pressure, etc. 

Below Tc and in zero field this Ising model has two degenerate magnetized phases corresponding to a majority 
of the spins in the positive or the negative direction. A weak applied field breaks the degeneracy, and the phase 
with the spins aligned (anti-aligned) with the field is stable (metastable) . If the field varies periodically in time, the 
system is driven back and forth across a first-order phase transition, and the two phases alternate between being 
momentarily stable and metastable. As a result, m{t) lags behind H{t), and hysteresis occurs. In the regime of small 
system size, weak applied field, and temperature well below Tc considered here, the system switches abruptly and 
stochastically between the two magnetized phases. A difference between two-dimensional, locally coupled bistable 
systems, such as this Ising model, and the one-dimensional arrays studied in most of the stochastic-resonance studies 
cited above |3^- |3g| , ^ , is that locally coupled one-dimensional systems have no ordered phase at nonzero temperature 
or noise intensity. The apparent long-range order in those studies is therefore a finite-size effect. However, the average 
equilibrium domain size grows exponentially with decreasing temperature p9| , ^^ . For chains much shorter than this 
size, the absence of true long-range order should not be qualitatively significant for the hysteretic behavior. 

The metastable phase in Ising models exposed to a static field H decays by different mechanisms, depending on the 
magnitude of H, the system size L, and the temperature T |Q. Two distinct regimes are separated by a crossover field 
called the dynamic spinodal, Hdsp{T, L). These two decay regimes can be distinguished by the statistical properties of 
the lifetime of the metastable phase. The lifetime is defined as T=t{m = 0), the first-passage time to a magnetization 
of zero, following an instantaneous field reversal from H to —H. For \H\ ^ H^sp, the mean of the lifetime, (r), 
is much greater than its standard deviation, ar- Therefore, this field region is termed the "deterministic regime." 
For \H\ <C i?DSP, (t) ~ and this field region is therefore termed the "stochastic regime." Both the deterministic 
and stochastic regimes are further subdivided according to the modes by which the metastable phase decays. The 
deterministic regime is split into the multi-droplet (MD) and strong-field (SF) regions for the low and high fields in 
this regime respectively. For a given system size the stochastic regime is also divided into the coexistence (CE) and 
single-droplet (SD) region for the low and high fields in this regime, respectively. Detailed discussions of these different 



decay modes are found in Refs. |45-47|. At sufficiently low T that the single-phase correlation lengths are microscopic, 
the different decay regimes can be distinguished by the interplay among four length scales: the lattice spacing a, the 
system size L, the radius of a critical droplet Rc, and the average distance between supercritical droplets Rq. The 
latter two lengths increase with decreasing field strength: Rc oc 1/\H\ and Rq oc exp[const./|i?|''^^], where d is the 
spatial dimensionality. Here we consider specifically decay in the SD region, which is characterized by 

a<i?c<i<i?o. (1.1) 

In this regime, the decay of the metastable phase proceeds by random homogeneous nucleation of a single critical 
droplet of the stable phase, which then quickly grows to take over the system. We have previously proposed [^-pT 46 



that this decay mechanism may apply to, e.g., barium ferrite particles in the 50-70 nm diameter range [p5[ . The 
crossover to the MD region corresponds to Rq ^ L. As a result, the dynamic spinodal depends asymptotically on 
L as Hr,sp{T,L) - (lnL)-i/('^-i). In the SD region the critical droplet is much smaller than the system itself, and 
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the crossover to the CE region is marked by Rc L. The corresponding crossover field, called the thermodynamic 
spinodal, therefore depends on L as Hthsp{T,L) ^ . In recent exploratory studies we have shown that 

the response of a kinetic Ising model to an oscillating field is qualitatively different for the MD and SD regions. The 
details of the response in the MD region will be described in forthcoming papers 

Theoretical studies of hysteresis have been performed for several models, using a variety of methods. These include 
various studies of models with a single degree of freedom, equivalent to mean-field treatments of the Ising model [ pi] 33 1 , 
Monte Carlo (MC) simulations of the spin-1/2 Ising model p^j27| , |5^ -|60t , and several 0{N) type models ||7|p5|-^p7|, 31|. 
These studies were performed with variations in the details of the simulations and in the model parameters. Most of 
them indicate that the average hysteresis-loop area, {A)—~{§ m{H) dH) , appears to display power-law dependences on 
the frequency and amplitude of H[t). However, there is no universal agreement on the values of the exponents, either 
experimentally or theoretically. For the Ising model, nucleation effects that would lead to a logarithmic frequency 
dependence have been proposed 26 . A mean-field model exhibits a dynamic phase transition in which the mean 
period-averaged magnetization, {Q) = (ijj /2'k){§ m{t) dt), changes from (Q) 7^ to { Q)—0 l32|. Such a dynamic phase 
transition has been suggested from MC simulations of a kinetic Ising model as well ]5 ^571 , |59| , |60| , |63t . 

The work presented in this paper differs from most past theoretical and numerical studies of hysteresis in two 
important ways. First, mean-field models do not take into account thermal noise and spatial variations in the order 
parameter, thus ignoring fiuctuations which may be important in real materials. Second, most previous investigations 
of hysteresis in Ising models have considered the frequency and amplitude dependence of quantities such as Q and A, 
without considering the manner in which the metastable phase decays. 

Considering the nucleation-based single-droplet decay mechanism, we find that the average hysteresis-loop area 
exhibits an extremely slow crossover to a logarithmic decay with frequency and amplitude in the asymptotic low- 
frequency limit. This crossover is sufficiently slow that the behavior can easily be misinterpreted as a power law 
over several orders of magnitude in frequency. We also show that the average loop area and the residence-time 
distributions for the system magnetization exhibit evidence of stochastic resonance, and we provide a connection 
between the characteristic decay time of the residence-time distributions and the power spectral densities of the 
magnetization time series. We find no evidence of a dynamic phase transition in the SD region. 

The rest of this paper is organized as follows. Section ^ supplies background information on the simulation of 
the kinetic Ising model. In Sec. [II some general properties of the time-series data are discussed. In Sec. [V the 
probability that the system magnetization does not switch sign during a period of the field, PnotC"^), is derived. This 
derivation is central to theoretical calculations throughout this paper. Section |^ presents theoretical calculations 
and MC simulation data for the residence- time distributions (RTDs). Also, we define and calculate the characteristic 
time of the RTDs and show its relevance to the low-frequency behavior of the power spectral densities (PSDs) of the 
time series, which are analyzed in Sec. VI. Section VII discusses the hysteresis-loop area , the correlation between 
the magnetization and the field, and the period-averaged magnetization. Finally, Sec. VIII contains a summary and 
conclusions. 



II. MODEL 



The model used in this study is a kinetic, nearest-neighbor Ising ferromagnet on a hypercubic lattice with periodic 
boundary conditions. The Hamiltonian is given by 

H = -J^S,5j (2.1) 

(ij) i 

where Si is the state of the zth spin and can have the values Si=±l, runs over all nearest-neighbor pairs, and 

runs over all N=L'^ lattice sites. The order parameter is the time-dependent magnetization per site, 
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1=1 



The dynamic used is the Glauber |Q single-spin-fiip Monte Carlo algorithm with updates at randomly chosen 
sites. The time unit is one Monte Carlo step per spin (MCSS). The system is put in contact with a heat bath at 
temperature T, and each attempted spin flip from Si to —Si is accepted with probability Q 

^Arl \ exp(-/3Ai;,) 
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Here AEi is the change in the energy of the system that would result if the spin flip were accepted, and (3 = l/k-oT 
where fee is Boltzmann's constant. It has been shown in the weak-coupling limit that the stochastic Glauber dynamic 
can be derived from a quantum-mechanical Hamiltonian in contact with a thermal heat bath modeled as a collection 
of quasi- free Fermi fields in thermal equilibrium [ |65| . 

In this paper, all numerical calculations are performed for (i=2, L=64 and T—Q.^Tc- This value of T is sufficiently 
far away from the critical temperature so that the thermal correlation length is small compared to the critical droplet 
radius and the size of the system. The system is subject to either an oscillating field, H{t) = — iJo sin(tjt), or to a 
constant field of magnitude -ffo- 

As discussed in Sec. |, the decay of the metastable phase in the presence of an external field H proceeds by nucleation 
of droplets of the stable phase . Figure || shows the metastable and stable phases as local minima in the free energy, 

F{m,H,T) ^ F{m,0,T)-mL^H, (2.4) 

of a nearest-neighbor Ising model on a 64 x 64 lattice at T = O.STc ^'^^ H=0 there are two degenerate 

equilibrium phases of magnetization ±meq(T'), separated by a free-energy barrier of height proportional to L'^^^ . For 
H = Hq = O.IJ the value of m near -1-1, where F has its global minimum, is the stable magnetization. The local 
minimum near m——l represents the metastable phase. The convex parts of the barrier represent a single spherical 
droplet of one phase embedded in the other. The droplet is a collective excitation p3| through which the switching 
proceeds, and the critical droplet is the droplet configuration corresponding to the local maximum of at a given 
value of H . For H < the stable and metastable phases are reversed. 

The average number of droplets of the stable phase that are formed per unit time and volume is given by the field 
and temperature dependent nucleation rate, 

/(i/(t),T)«i?(r)|i?(i)|^exp 

The notation follows that of Ref. |ll9|], where B{T) is a non-universal temperature dependent prefactor, and K and 
So(T) are known from field theory pSj-fTOt and simulations |Q and are listed in Table |. The quantity So(r) is the 
field-independent part of the free-energy cost of a critical droplet, divided by ksT. The external field, H{t), is the 
only quantity through which I (H(t),T) depends on time in this adiabatic approximation. 

Several quantities, whose values do not depend on the frequency of the applied field, are required as input for 
the theoretical calculations in the following sections. These quantities, which include the average lifetime {t{Hq)) « 
[L'^I{Ho,T)]~^ , are listed in Table |. They are determined through what we refer to as "field-reversal simulations." 
In these simulations the system initially has all spins up or positive. It is then subjected to a static external field of 
magnitude Hq with a sign opposite the system magnetization. This instantaneous field quench prepares the system in 
a metastable state, and the decay of the metastable phase proceeds by the mechanisms outlined in the introduction. 



So(T) 



\H{t)\ 



d-i 



(2.5) 



III. TIME-SERIES DATA 



In the simulations presented here, a sinusoidal field is applied to the system. Its amplitude, Hq = O.IJ < i?DSPi 
is chosen such that in field- reversal simulations the system is clearly in the SD region for a field of magnitude Hq. 
The dynamic spinodal field is approximated by i?DSP ~ ^1/2 j where H1/2 is the value of H (for given L and T) for 
which the relative standard deviation of the lifetime, r=ar/ (t), is 1/2 ES]. This value of Husp is given in Table |. 
It approximately equals the field for which the local minimum in Fig. |^ disappears 66 1. 



To obtain the raw time-series data, an Ising system was initially prepared with either a random arrangement of 
up and down spins with mit = 0) « 0, or with a uniform arrangement with all spins up. Then the sinusoidal field, 
H{t)——HQsm{Lot), was applied and changed every attempted spin flip, allowing for a smooth variation of the flcld. 
The time series did not appear to depend on the initial conditions after a few periods. The simulations were performed 
with several values of the driving frequency lo. For each frequency, we recorded the time-dependent magnetization 
m{t) for approximately 16.9 x 10^ MCSS. Each of these raw time-series data files store the values of t, H{t), and m{t) 
in increments of 1 MCSS. Each file takes up about 800 megabytes and took about 9 days to run, using a single node 
of an IBM sp2. These are by far the most extensive MC simulations of hysteresis in Ising systems to date. 

It is useful to think of hysteresis as a competition between two time scales: the average lifetime of the metastable 
state following an instantaneous field reversal from Hq to —Hq, {t{Hq)), and the period of the external forcing field, 
2tt/oj. Therefore, we specify the ratio R of the period to the average lifetime, 

p_ (W^) .o^^ 
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One may think of i? (1/i?) as a scaled period (scaled frequency). 

We note that {t{Ho)) is the "shortest of the long time scales" in the present system. From Fig. |l| we observe 
that whereas the free-energy barrier represented by the critical droplet for \H\ — Hq is on the order of ksT, the 
barrier between the degenerate phases at = is on the order of bOksT. For this temperature and system size, 
the timescale for spontaneous fluctuations between the phases in the absence of an applied field, (t(0)), is therefore 
essentially infinite. Conversely, the nuclcation of the critical droplet necessary to leave the metastable phase is 
entirely driven by the thermal fluctuations, even when the field has its maximum strength Hq. Switching in this 
system therefore truly depends on the joint action of the random thermal noise and the deterministic oscillating field. 

Figure ^ shows short initial segments of the magnetization time series in the SD region for three different values of i?. 
In all three cases, mit) fluctuates near one of the two degenerate values of the spontaneous zero-field magnetization, 
punctuated by rapid transitions between these two values, that are completed during a single half-period of the 
applied field. The rapid switching of m{t) is evidence of the nuclcation of a single critical droplet that reverses the 
sign of the magnetization. The magnetization 'plateaus' are due to the failure of any critical droplet of the stable 
phase to nucleate. The fluctuations in m(t) on these plateaus indicate appearance and disappearance of sub-critical 
droplets. For low driving frequencies the magnetization switches twice during almost every field cycle, whereas for 
high frequencies switching occurs only occasionally. The 'spikes' in m{t) seen for i?=2.5 occur when a single droplet 
nucleates, but does not have time to grow and switch the system magnetization before the field becomes unfavorable 
and the droplet rapidly collapses. The number of field cycles shown in Fig. ^ is small compared to the total number 
of cycles in an entire time series. 



IV. PROBABILITY OF NOT SWITCHING DURING A PERIOD 



The probability that the system does not switch during a full period of the field, Pnot(w), is central to the theo- 
retical understandingof hysteresis in the SD region. It occurs most directly in the calculation of the residence-time 
distributions in Sec. |V[ In addition, elements in the derivation of Pnot(^) are fundamental for describing most of the 
observed quantities in this study. 



As mentioned in Sec. Ill, the system exhibits abrupt switches during which the average magnetization changes 
between values near ±mcq. As seen in Fig. ^, these events occur quickly compared to the period of the external field, 
and to a first approximation the time it takes the droplet to grow to fill the system (the growth time) is negligible. A 
more realistic treatment takes the finite growth time into account as a lag time between the nuclcation of a critical 
droplet and the time at which the system switches. 

The first part of the derivation is presented without the effects of the growth time. This is done for simplicity, as 
well as to emphasize the role of the growth time as a correction to the basic picture of a variable-rate Poisson process. 
First we derive the expression for the cumulative probability that a switching event has occurred by time t, F{t), in 
terms of the time-dependent rate of an instantaneous decay process, p{t). [This cumulative probability should not be 



confused with the free energy F(m, H, T) of Eq. (2.4).] It is convenient to introduce F{t) = 1 — F{t), the probability 
that a switching event has not occurred by time t. Standard theory of variable-rate Markov processes |Q leads to a 
difference equation for F{t), 

F{t + At) = F{t) X {probability an event has not occurred in the interval [t, t + At]} 

^F{t)[l-pit)At] , (4.1) 



which in the limit At gives 



dF{t)/dt^ -p{t)F{t). (4.2) 



The growth time, ig(i), is introduced into the derivation at the level of Eq. (4.1). It is defined as the time between 
the nuclcation of a critical droplet and the time when the volume of this droplet becomes approximately half the system 
volume. The dependence of the growth time on t is a consequence of the time dependence of the interface growth 
velocity, which is approximately proportional to H{t). For suitably long time scales, the growth of a supercritical 
droplet is a deterministic process. Another quantity in this derivation is the time at which a droplet nucleates, tn{t). 
If a switching event occurs at time t, then tn{t)=t — tg{t). Where clarity is not sacrificed, we do not show the explicit 
t dependence of tg or t^- For a switching event to occur in any particular period of the external field, a critical droplet 
must not only nucleate, but must do so early enough so there is sufficient time for it to grow to the volume of the 
system. Therefore, the difference equation for F{t) is modified to read 
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F{t + At) = F{t) X {probability a switch has not occurred within [t, t + At]} 

— F{t) X {probabihty a droplet has not nucleated within [tn,tu + At^]} 
= F{t) [1 - p(t„)At„] , 



l^pit^)^At 



We can express this result in terms of the growth time tg and its derivative, using 

dtn ^ dtg 
lit " ~ Ht ' 



Substituting into Eq. (4.3) and letting At ^ gives 

dFit) 



dt 



Pit ~ <g) 



1 



dtg 
dt 



F{t). 



Integrating Eq. (4^) gives the cumulative distribution, 

■ dhit')^ 



F{t) = 1 - exp 



P[t'-kit')]{l-^)dt' 



Differentiation gives the probability density function (pdf) for switching events at time t, 



Pit) = ^ = p[t-h{t)] 



dt 



exp 



p[t'~tg{t')](l 



dt' 



dt' 



(4.3) 



(4.4) 



(4.5) 



(4.6) 



(4.7) 



For tg—0, Eq. ( ^77| ) is equivalent to Eq. (9) of |7^. The growth time tg is obtained from the expression for the 
time-dependent volume of a supercritical droplet 



v{t,t^)^n 



v{t')dt' 



tn 



(4.8) 



where t > t^- Here v{t) is the droplet interface velocity and O is defined such that the volume of an equilibrium droplet 
of radius R is flR'^ ||7l|. Using the Lifshitz-Allen-Cahn approximation [[75[-|78|, the interface velocity is v{t) w i^\H{t)\. 
The proportionality constant v depends on the details of the dynamics. Here we use values for the Glauber dynamics, 
obtained from field- reversal simulations by Ramos et al. fl% . The values of the constants used in our calculations are 
listed in Table |. For d^2, Eq. ( |4.8D for the growth time becomes. 



v{t,t-tg) = — = n 



For a static field of strength ffo, the growth time is 



t. = 



1 



vH^ sinLot'dt' 



|coscj[t — tg] — cosoj^l 



L 



Substituting this expression into Eq. (4.9) gives 

cosa;[t — tg] = tgLu + cosut. 

Solving for tg such that tg < t gives 

t — cos^^[cosujt + igu;]} tQ<t<TT/uj 



kit) 



otherwise 



(4.9) 



(4.10) 



(4.11) 



(4.12) 



where 
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^0 = — COS "'"[1 — tgCj] 



(4.13) 



The time to is the first time during a period for which the probabihty of switching is non-zero. 
m{t=0) ~ — 1, then the probabihty density P{t) that a switching event takes place at time t is 



If H{t^O)^0 and 



Pit) 








dt 



exp(-/>[i'-tg(t')] 



rftg(t') 

dt' 



dt 



< < < to 
to <t < tt/w 
tt/uj <t< 27r/w, 



(4.14) 



where the ranges for t ensure P{t) is non-zero only when t > tg and the signs of ■m{t) and H{t) are not equal. Higher- 
order corrections, including the probability that a second droplet nucleates during tg, were found to be numerically 
insignificant. The main approximation used here to obtain tg{t) lies in ignoring the slower growth of droplets only 
slightly larger than the critical radius. This has been shown to be permissible for adiabatically slow- forcing models in 
which large droplets grow exponentially in time fz^ . In the present case, however, we simply consider it a convenient 
approximation, whose accuracy is ultimately confirmed by our numerical simulations. 

In the SD region, the average lifetime in a field-reversal simulation should be dominated by nucleation. Therefore, 
the total nucleation rate in a static field Hq can be expressed as 



= L''liHo,T). 



(4.15a) 
(4.15b) 



The nucleation rate in a static field of strength Ho should be equal to the nucleation rate in a sinusoidal field of 
amplitude Ho at the maximum of the field, po=p{t = Tr/2uj). The ratio of these two decay rates is then 



Po 



I{H{t),T) 



(4.16) 



Substituting the form of the nucleation rate, Eq. ( p. 51) , into the expression above allows p{t) to be recast in a form 
which does not explicitly contain the non-universal prefactor B{T): 



p{t) = Pol sin(a;t)|''^' exp 



So(T) 



1 



sin(wt)|''- 



- 1 



(4.17) 



This expre ssion holds when ■m{t) and H{t) have opposite signs, while p{t)~0 when they have the same sign. Usin, 
Eq. ( 4.15a ) for the maximum decay rate gives po=(6.62± 0.07) x lO^"' MCSS^^, using quantities listed in Table 
Figure |^ shows P{t)/uo vs. uot for five different frequencies of the external field. The inset in Fig. |3| shows p{t) vs. ut 
The nucleation rate achieves its maximum value, poj at a phase of ut—'K/i^ independent of w. However, the location 
and width of the maximum for P{t) depend strongly on uj. This behavior results from the combined field dependence 
of the nucleation rate and the interface growth velocity. For R > 100, P{t) narrows and the location of its maximum 
shifts to lower phase values as the switching begins to occur before the maximum in p{t). As R is decreased below 
1.5, ujto TT, and the area under the curve for P{t)/uj goes to zero. Therefor e, th e condition Wmaxto=7'' gives the 
maximum frequency for which single-droplet switching is possible. Using Eq. (4.13) and converting the result to a 
bound on 1/i? gives (l/i?)max=l-19 ± 0.03. For higher frequencies, switching events are very rare, and if they occur 
at all, they do so through a multidroplet mechanism. The probability of not switching during an entire period is 
obtained by integrating the probability density P{t), 



Pnot{^) = l- P{t')dt' 

Jo 



(4.18) 



The frequency dependence of Pnot(w) is the aspect of this quantity most important for comparisons with our MC 
data. However, Pnot(w) also depends on other parameters through the nucleation rate. 



V. RESIDENCE-TIME ANALYSIS 



As mentioned in Sec. Ill, the magnetization exhibits abrupt switches between values near ±meq. In contrast, the 
times between the magnetization reversals are comparable to, or greater than, the period of the applied field. As a 
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first approximation, one may therefore consider these switching events as occurring in a discrete two-state system. 
For an Ising system undergoing a field-reversal experiment ||45|] , the lifetime in the SD region is stochastic and is 
well described by droplet theory. For an oscillating field, the analogous quantity is the time between reversals of the 
magnetization, called the residence tim e. The probability density for the residence times is called the residence-time 
distribution (RTD) In Sec. ^ we construct an alytic al expressions for the RTDs and compare these with the 

RTDs obtained from our simulated time series. In Sec. VB we calculate the area of the peaks in the RTDs, or the 
RTD peak strengths, and compare our theoretical results for the peak strengths with MC data. Finally, we show that 
our data for the RTD peak strengths provide evidence of stochastic resonance in the model. 



A. Residence-time distributions 



We define the residence time. A, as the time between consecutive magnetization reversals, and denote its probability 
density as n(A). The details of our theoretical deviation of n(A) are given in Appendix ^ The results of the 
theoretical calculation for the residence-time distributions, which contain no adjustable parameters, are shown as 
solid curves in Fig. ^(a-c) for different values of R. 

Next, we give a description of the MC analysis for the RTDs. (The results of this analysis are shown as solid dots 
in Fig. ^a-c) for different values of R.) When measuring an RTD one must ignore "false crossing events." In these 
events, the magnetization crosses zero and re-crosses zero again within a short time without having reached a value 
near the stable magnetization. There appear to be two reasons for these re-crossing events. First, when m(t) w the 
magnetization can re-cross zero many times due to thermal fluctuations. Second, during some of the periods there 
will be a "spike" in the magnetization when a supercritical droplet nucleates but does not have time to completely 
take over the system before the applied field changes sign. For this reason, a cutoff is employed, and a switching event 
is recorded only when m{t) reaches some cutoff value ±mcut- To quantify the meaning of the cutoff, define t^' and t~ 
as the times at which m{tf)=+Tncut and TO(i^)=— mcut, respectively. The residence time A^ is given by 

^ r 1 - t- when m{t) w -1 for <t< t+^^ , , 

' \ tr^^ - tj when m(t) w +1 for t+ < t < t^^^ ' ^ ' ' 

We used TOcut=0-25. For each frequency, the residence times are measured over an entire time series. The size of the 
bins in an RTD is set by dividing the maximum observed residence time by the number of bins. Both the maximum 
residence time for a given time series and the number of bins are different for different frequencies. Hence, the size 
of the bins is different for each of the graphs in Fig. |4[ Scaling the residence times by the period of the external field 
centers the peaks in the RTD about every odd half-integer. In the low-frequency limit, the system spends enough 
time in an unfavorable field during every half-period to allow the magnetization to switch. In this limit, the RTD 
would contain a single peak centered around 1/2. As the frequency of the field is increased, there should be more 
periods during which the magnetization does not switch at all, indicated in the RTDs by an increase of the size of the 
peaks centered on 3/2, 5/2, etc. The RTD data from MC simulations are shown in Fig. ^ as solid points together 
with the theoretical curves. 

The smoothest MC results, and the best agreement between theory and simulation occur for R ^ 10. The agreement 
is quite good, considering that the theoretical calculation contains no free parameters. All of the constants in the 
formulas for the nucleation rate and the interface growth velocity come from theoretical considerations or from field- 
reversal simulations. The agreement is poor only for the highest frequencies, corresponding to i? ^ 2.5. For these 
values of R the MC data are suspect. First, for these high frequencies of the external field, the sizes of the peaks 
for large residence times are significant. For all of the frequencies shown, the data sets have approximately the same 
total number of switching events. Therefore, a smaller number of events is contained in each bin of the RTDs for the 
higher frequencies. Second, in spite of the cutoff there are two peaks in the RTDs for R=2 through 5 in the interval 
< ujA/2tt < 1. Of these two peaks, the peak at the shorter residence time comes from "spikes" in the magnetization 
which are large enough to extend past the cutoff value, but still do not switch the system completely within a field 
period. These "spikes" in the time series redistribute weight from the higher-order peaks of the RTD into the spurious 
peak at short residence times. This affects the measurement of the peak strengths from the MC data as well, as 
discussed in the next section. 



B. RTD peak strengths 



The frequency dependence of the strengths of the peaks in the RTD is another quantity which describes the nature 
of the magnetization reversal. The strength of the j-th peak in an RTD is given by 
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Sj{uj) = / n(A)dA. (5.2) 

•^0-1)^ 

The peak strengths obtained from the MC time series are shown as sohd dots in Fig. |5[ The statistical errors are 
estimated by error-propagation analysis as Sj{l — Sj)/N, where N is the total number of switching events in a 
simulation run. For almost all of the data points the error bars are smaller than the symbol size. The solid lines in 
Fig. H are theoretical results. The strength of the first peak, 5'i(w), is simply the probability that m{t) switches sign 
within the first period after the last field reversal, 

Sliiv) ^ I - Pnotiuj) . (5.3) 

Therefore, the strength of the j-th peak is 

Sj{u;) = Pnotio^y-' [1 - Pnot(^)] ■ (5.4) 



The values for Pnot('^) used in this calculation were obtained by numerically integrating Eq. (4.18) for several values 
of Lo. This parameter- free numerical evaluation of the theoretical peak strengths is in good agreement with the MC 
data. In Fig. || one can see this agreement, especially for the strength in the first peak, 6*1(0;), for all but the highest- 
frequency data point at l/i?=0.5. However, the MC data slightly overestimate 5*1(0;) even for low frequencies. This 
is due to the redistribution of strength from the higher-order peaks into the first peak due to "spikes" in the time 



series which extend past the cutoff, as mentioned in Sec. VA. Hence, the peak strengths for the higher-order peaks 



are systematically underestimated by the MC data, particularly for 6*2(0;). The agreement between the theoretical 
curve and the data is not quite as good for 62(0;), 6*3(0;), and 64(0;) at higher frequencies. However this is expected, 
due to the poorer statistics for these higher-order peaks. 

Analysis of the RTDs for two-state systems has been used to detect stochastic resonance (SR) | |30| , |80t . Gammaitoni 
et al. studied the switching behavior and residence times of an analog circuit, which served as a model of a bistable 
system driven by random noise and a sinusoidal external forcing . They found for their model that SR is manifest 
in the fact that each of the peak strengths in the RTDs has a maximum for a given frequency of the field. If the 
frequency Uj corresponds to the location of the maximum in Sj , then 

o;i < o;2 < o;3 < . . . . (5-5) 

This approach gives an alternative definition for SR, different from the original definition as the noise intensity for 
which the signal-to-noise ratio (SNR) exhibits a maximum pO| ]. Since the SNR does not exhibit a maximum with 
respect to the frequency of the forcing, the definition suggested by Gammaitoni et al. facilitates the understanding 
of SR as the tuning of one timescale (inverse frequency) to another (average lifetime of the metastable phase) , more 
analogous to a "bona fide" resonance. 

The case studied in Ref. |^l[] was one of a very weak oscillating field, such that the noise-driven switching rate in 
zero field was of the same order of magnitude as the escape rate from the metastable state in maximum field. As 
a result, there was only one long timescale, and for sufficiently low frequencies most of the escapes were completely 
thermally driven. This would redistribute the number of escapes with residence times less than n/iv into a peak at 
much shorter times, corresponding to the thermal escape rate of the unforced bistable system and giving a nonzero 
value of o;i . 

For the system studied here the thermal switching rate in zero field is virtually zero, as pointed out in the discussion 



following Eq. (3.1). As a result, 6*1 increases monotonically towards unity as the driving frequency is lowered. For 
61 to display a maximum, the weight in the RTD's must shift towards times much shorter than the period of H(t). 
These residence times would correspond to events in which a single thermal fiuctuation switches the system. One 
can increase the thermal switching rate in zero field (or decrease the thermal relaxation time in zero field (t(0))) by 
increasing T. (Recent measurements of telegraph noise in nanoscale ferromagnetic particles are able to show this quite 
clearly [p2[.) But increasing the temperature can move the system into an entirely different decay regime (either the 
MD or SF regime, depending on the value of T) where the stochastic nature of the response disappears. However, 
even thou gh 6 *1 does not display a maximum, the higher-order peak strengths do have maxima at o;2 < o;3 < 



From Eq. (5.4) the theoretical positions of these maxima are seen to be given by the condition Pnot(wj) = {j ~ 
which yields o;i = 0. Through i-'„ot(o;) they are determined by the competition between the period of the deterministic 
forcing and a stochastic timescale, which is the minimum metastable lifetime {t{Ho)). 

If we were to reduce L to push the simulations into the coexistence (CE) region, then (t(0)), which depends 
exponentially on L"^^^ through the barrier in the free energy F{m,0,T), would approach the escape time {t{Ho)), 



which by Eq. (4.151) is inversely proportional to L'^. As a result, 5*i should be observed to decrease at very low 



frequencies for sufficiently small systems. In this regime, the critical droplet volume would be on the order of half 
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the system volume. This effect was recently observed by Lindner et at. in simulations of a one-dimensional chain of 
bistable elements driven at a constant frequency and subject to noise of variable intensity [|3 6(1371 ] . 

Our results lead us to make the observation that a maximum in Si is not necessary for the response of the system 
to the oscillating field to be characterized as "resonant." Rather, m{t) is essentially synchronized with H{t) in the 
whole frequency range where 5*1 is close to unity. The upper limit of this range is proportional to L02 and therefore 
determined by {t{Ho)), whereas the lower li mit w ould be determined by the (in our case unobservably long) (r(0)). 
This theme will be discussed further in Sec. VII A. 



C. Characteristic time of the RTDs 



For each of the RTDs shown in Fig. ^, the size of the peaks decreases for large residence times. The rate of this 
decrease can be quantified by measuring how the peak strengths decrease with increasing peak number. We define 
the characteristic time, r]{ui), for the RTDs (in units of the field period, 2tt/ui) as 



1 



In Sj (w) — In Sj+i{uj) 



(5.6) 



The value of ?/(aj) for any frequency of the field is measured from the MC data by plotting InS'j vs. j. The slope of 
a best- fit line through this data gives —1/77(0;). The inverse characteristic times calculated from the MC data by a 
weighted least-squares fit are shown in Fig. ^ The statistical uncertainty in the estimates for the characteristic time 
becomes large for the lowest and highest frequencies. For the very lowest frequency shown, l/i?=0.05, few switching 
events contribute to any of the peaks in the RTD, other than the first peak. Hence there are poor statistics in the 
InS*-,- data for j > 1, resulting in a large uncertainty in the fitted slope. For high frequencies the RTDs contain many 
peaks. Since the total number of switching events in each time series is approximately equal, the statistics for each 
peak is poor. In addition to this purely statistical error, there is also a systematic effect due to the cutoff used to 
measure the residence times. Namely, the RTDs display two peaks for the interval < ujA/2tt < 1. This extra weight 
in Si{uj) introduces additional systematic error which tends to raise the estimates of 1 1/77(0;) | for higher frequencies. 

The theoretical calculation of the characteristic time starts by substituting Eq. (5.4) for Sj{uj) into Eq. (p!q), which 
gives 



77(0;) = 



-1 



InPnot(t^) 



(5.7) 



The calculation of 77(0;) is now trivial because Pnot('^) has already been evaluated numerically in Sec. VB, The solid 
curve in Fig. ^ shows the theoretical results from the full numerical calculation of Pnot(w). The theoretical result 
diverges at a value of 1/i? « 1.19, the maximum value of the frequency for which single-droplet switching is possible. 
At this value of 1/R, Pnot=l, and the characteristic time diverges. For frequency values 1/R > 1.19 mechanisms other 
than single-droplet decay might be possible, such as several droplets nucleating simultaneously. Given the length of 
our simulations, it is unlikely such behavior would be observed. Indeed, the statistical and systematic errors discussed 
above preclude accurate measurement of the characteristic time from the MC data for 1/R> 0.4. 

For low frequencies, both the theoretical and the MC results for the characteristi c ti me appear constant. This 
motivates the search for an approximate, analytic expression. We start by casting Eq. (5.7) in terms of the cumulative 
probability distribution. 



Using the low- frequency approximation tg=0 gives 



(5.8) 



r){uj) 



2/- pit')dt' 



(5.9) 



where the upper limit in the in tegra l is rewritten because the nucleation rate p{t) is symmetric about i=7r/2o;. 
Substituting it=sin(o>f) into Eq. (4.17) for p{t) allows one to write the integral above as 



77(0;) 



1 

2po- 



: exp 



Ho \u 



du 



(5.10) 
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The final result for the characteristic time is given in units of (r) by multiplying both sides of the equation above by 
R and simplifying, 



Rf] 



Po ' 



: exp 



Ha \u 



du 



(5.11) 



where all the quantities are known except for the integral, which can be calcula ted numerically. This theoretical result 
is shown in Fig. ^ as a horizontal dashed line. The definition of ri{u:) in Eq. (x6) admits to an interpretation of the 
characteristic time as an exponential decay constant for Sj(uj) as a function of j. Therefore, ijiuj) is a convenient 
average measure of an RTD, with long characteristic times corresponding to large peaks at long residence times. 
The exponential nature of the decay of the RTD with time should give rise to a Lorentzian component in the power 
spectral density with half width fi/2 ~ [2TrRr](r)]^^ w 1.4 x 10~^MCSS~^, obtained from a numerical evaluation of 



Eq. (5.11). This low- frequ ency component should be visible in the spectrum when fi/2 < uj/2n. 

The result in Eq. (5.11) was obtained in a low-frequency approximation by setting tg—0. Further approximation 
can be made in the evaluation of Eq. (5.11). The main contribution to the integral in the denominator occurs for 
M « 1. Expanding each factor in the integrand in e=l — u for small e and substituting x'^—'E.Q(T)e/ Hq gives the 
Gaussian approximation. 



Rt] 




(5.12) 



which gives [2TTRr] (t)]"^ =1.9 X 10-5MCSS"\ and was obtained by inserting quantities from Table |. The agre emen t 
between this value and the data is not as good as that obtained directly from the numerical evaluation of Eq. (5.11). 
However, the expansion in small e improves for large So(T) / Hq. So this analytic formula for r]{uj) would be appropriate 
for low frequencies and small amplitudes of the external field, i.e. field amplitudes which place the system even deeper 
into the SD region. 



VI. POWER SPECTRAL DENSITIES 

A standard method used to characterize a time series is to calculate its power spectral density (PSD). Figure 
shows the PSDs of the raw data, short segments of which are shown in Fig. ||. These PSDs are calculated with a 
standard FFT algorithm implemented using a Welch window |Q. To reduce the variance in the PSDs, each time 
series is split into several segments. The data is then overlapped in such a way as to obtain an "average" over all the 
segments for each frequency bin. The details of this method, which we refer to as "smoothing," are given in Ref. [^ . 
Depending on the number of segments into which the original time series is split, there is a trade-off between frequency 
resolution and variance reduction. If the time series is split into many segments, the variance per bin is decreased at 
the expense of lower frequency resolution. Conversely, by partitioning the data into fewer long segments one obtains 
higher frequency resolution and a smaller low-frequency cutoff, at the expense of a larger variance per bin. 

The PSDs for different driving frequencies are shown in Fig. ^ The spectra in the main part of the plot have been 
shifted in the vertical direction by arbitrary offsets. The spectra in the inset of Fig. |^ are plotted with no offset. 
Different amounts of smoothing have been used for the low- and high-frequency regimes. In the low-frequency regime, 
less smoothing is used. This increases the frequency resolution, enabling one to see sharper peaks in the spectra and 
sample more of the low-frequency response. In the high-frequency regime, more smoothing is used. Since there is less 
power in the PSDs for the higher frequencies, a reduction in the variance is a higher priority than frequency resolution. 
The fourth spectrum shown in Fig. 0, labeled "background," corresponds to thermal fluctuations in a single-phase 
system. To obtain this spectrum, a simulation was performed on a system with the same size, temperature, and for 
the same number of MCSS as the other spectra, in a static field of Hq / \/2. The magnetization quickly relaxes to 
equilibrium. The equilibrium fluctuations are purely thermal, and their time correlations are exponential with a short 
correlation time of only a few MCSS. The PSD should then have the functional form of a Lorentzian. When plotted 
on a log-log scale a Lorentzian appears flat at low frequencies, then crosses over into a linear curve with a slope of 
—2. The tail on the observed noise background does not appear to have slopc=— 2. However this is most likely due 



11 



to aliasing effects near the Nyquist frequency, ^n—t^ which is only about one order of magnitude larger than the 
inverse correlation time. 

To describe the PSD for each frequency, we identify four distinct regions: 1.) the thermal noise region, 2.) the 
peaks, 3.) the low-frequency region, and 4.) the intermediate region. The inset in Fig. ^ shows how the PSDs 
collapse onto the background spectrum in the thermal noise region for large frequencies. The timescale of the thermal 
fluctuations is much shorter than the shortest period of the external field. 

The most prominent features of the PSDs are the sharp peaks. For each driving frequency, the first peak in the 
spectrum is located at w, the frequency of the external field. The second peak is located at Slo, the third peak is 
located at Sw, and so on. These odd harmonic peaks arise because the shape of the time series strongly resembles a 
square wave (refer to Fig. |2). The power p„ contained in the nth component of the Fourier series for a pure square 
wave is p„ — 16[sin(7i7r/2)] /(n7r)^, which is non-zero only for odd n and decays as n~^. The zeros in p„ for even n 
survive as zeros in the PSD for a switching process in which the switching probability density P{t) reduces to a delta 
function [Q. The finite width of P{t) for the present system smoothes out such singularities in the PSD. However, a 
clear dip in the PSD at twice the driving frequency can be seen for R=10. This effect was also noticed by Zhou and 
Moss in the weak- noise regime for analog simulations of a bistable circuit |^ . 

The low-frequency region comprises the portion of each spectrum between the first peak and the lowest resolved 
frequency. The PSD in this region exhibits a strong dependence on the frequency of the external field. For _R=100, 
the average intensity in the low-frequency region is approximately constant and is weaker than the fundamental peak 
intensity by about three orders of magnitude. For R—2.5, the slope of the PSD in the low-frequency region is close to 
—2 over almost two orders of magnitude and contains components comparable in intensity to that of the fundamental 
peak. This significant amount of power at low frequencies is a consequence of the long residence times, i.e. those 
residence times longer than the period of the field. The Lorentzian half width predicted from the RTDs in Sec. |V C| , 
fi/2 ~ 1-4 X 10~^MCSS~^, is in good agreement with the PSD for i?=2.5. One can see from Fig. |7|that if the first peak 
for a PSD is located at a frequency smaller (larger) than approximately 1.4 x 10~^MCSS~"'^, there are small (large) low- 
frequency components. In an analogous fashion, the low-frequency behavior of the PSDs can be deduced by comparing 
the inverse characteristic time, [2iTRri (r)]~^, in Fig. ^ with the frequency, (t) / R (in units of MCSS~^), which is shown 
as a short-dashed line in the same figure. For low driving frequencies, (r) /R < [2nRr] (r)] ^ . Therefore, most of the 
residence times are less than a period of the external field, which corresponds to small low-frequency components in 
the PSDs. For high driving frequencies, (t) / R. > [27ri?77 (r)] , and large low-frequency components appear in the 
PSDs. 

The intermediate region is the portion of each spectrum between the highest-order visible peak and the thermal-noise 
region. This region is discussed last because it is best understood as a crossover from the peaks to the thermal-noise 
region. The structure of the odd-harmonic peaks is easiest to discern for i?=100, where the first five peaks are clearly 
visible. The higher-frequency harmonics cannot be seen for two reasons: the peak positions become very closely 
spaced on a logarithmic scale, and the intensity of the peaks becomes too small to be seen above the fluctuations in 
the PSD. So the intermediate region may be thought of as an envelope of the high-frequency, odd-harmonic peaks. 
A log-log plot of pn vs. odd n yields a line with a slope of —2. However, the portion of the PSD just to the left 
of the thermal-noise region is decreasing with a slope steeper than —2. This sharp falloff is a consequence of the 
flnite growth time of a critical droplet. A flnite growth time effectively "smoothes out" the sharp corners of the 
square-wave response of the magnetization, which introduces a cutoff for the highest-frequency components of p„ in 
the intermediate region. 



VII. HYSTERESIS LOOPS, CORRELATION, AND PERIOD-AVERAGED MAGNETIZATION 

To characterize the behavior of an entire time series we calculate the following integrals of the magnetization. 



A = -<l)m{H)dH (7.1) 



— i m{t) H{t) dt (7.2) 
27r / 

Q = — i m{t) dt , (7.3) 
2vr J 

where A is the hysteresis-loop area, B is the correlation between the external held and the system magnetization, and 
Q is the period-averaged magnetization. These quantities are calculated over each period in the entire time series. 
From the resulting "flltered" time series we construct histograms to obtain the probability densities of A, B, and Q 
for each separate frequency of the external fleld. 
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In Sec. VII A we compare our MC data and theoretical calculations of A and B. In particular, we comment on the 
low- frequency power-law scaling for A which has been put forth in numerous studies. At the end o f this section we 
show our results for B and identify A and B as components of a nonlinear response function. In Sec. VII B we discuss 
Q and the absence of a dynamic phase transition in the SD region for this system. 



A. Hysteresis loops and correlation 

The hysteresis-loop area, A, represents the energy dissipated during a single period of the applied field. It is therefore 
one of the most important physical quantities characterizing hysteretic systems, and it is frequently measured in 
experiments. Under the conditions of stochastic magnetization reversal, A and B are random variables with nontrivial 
statistical properties. Figure ^ shows the probability densities of A. For all values of R ^ 10, there is a sharp peak 
near zero. This peak is denoted as Aq and corresponds to the field cycles during which m{t) does not switch sign, but 
merely fluctuates near ±meq. The second peak, Ai, is located near A/{AHo)=0.5. It represents field cycles during 
which m{t) switches sign once. The third peak, A2, located near a loop area of A/(4i?o) = l-0, represents cases when 
■m(t) switches sign twice within the same period. When the period of the field increases, the weight in the peaks moves 
from low to high values of A. For H{t) with longer periods, the magnetization has a higher probability of switching 
once or twice during a single period, thus transferring more weight to the peaks Ai and A2. For very low frequencies, 
the magnetization almost always switches twice in each period, giving loops of type At the same time, switching 
occurs earlier in each half-period. This reduces the average loop area, giving rise to the type of pdf shown for i?=100 
in Fig. ||. Figure^ shows typical hysteresis loops with values of A corresponding to the three peaks, Aq, Ai, and 

The fraction of events contained in each of these peaks we denote as /o, fi, and /2. These are the probabilities that 
the magnetization switches zero, one, or two times in a full period of the field. They depend on the probability that 
the sign of m(i) is opposite that of H{t) immediately after H{t) has switched sign at the beginning of the period. We 
call this phase probability p-. In terms of p- and Pnot(w) we have: 

/o = Pnot(^) (7.4a) 

/l =P-P„ot(c^)(l - Pnot(^)) + (1 - i^not(^)) (7.4b) 

/2 =P-(l-F„ot(^))' ■ (7.4c) 

The phase whose probability is given by p- changes in a two-state Markov process described by a transition matrix 
fully determined by Pnot(ijj)- It is easy to show that the stationary value of p-=l/(l -l-Pnot(w)) After the phase 
distribution has reached its stationary state, we therefore get 

/o = Pnot(w) (7.5a) 
_ 2Pnot(u;)(l-Pnot(u;)) . 

(1 + Pnot('^)) 

(l-Pnot(^))" , . 

= JTTKM) • ^'-''^ 

These expressions satisfy the normalization condition /o + /i + /2 = 1- The fi are directly related to the RTD peak 
strengths, Sj, through their common dependence on Pnot- In particular, /o=Pnot=l ^ Si. 

It is less clear how to separate the peaks when measuring the fractions from the MC data. For example, with cutoffs 
at 0.333 and 0.666 an event is considered part of peak if < A/AHq < 0.333, part of peak Ai if 0.333 < A/AHq < 
0.666, and part of peak A2 if 0.666 < A/AHq < 1.0. Different cutoffs were tried and judged by how well they fit the 
theoretical results. No particular choice produced a fit that was qualitatively better than others. Comparisons of 



the fi calculated from theory and simulation are shown in Fig. 10, The data points are measured from the loop-area 
distributions, using cutoffs at 0.333 and 0.666. The solid curves are calculated from Eqs. (7^a-c), using the theoretical 
values for Pnot('^)- Considering that no free parameters are used in the theory, it agrees well with the simulation data. 

The average loop area for a specific frequency is a quantity often displayed in experimental and numerical studies of 
hysteresis. It can also be obtained in the present case, and we do so below. However, for stochastic hysteresis this is 
not a very useful quantity at higher frequencies, where the distributions are trimodal. The means of the distributions 
in Fig. H are shown in Fig. |l^. Note that the vertical bars are not error bars, but rather give the standard deviations 
of the distributions. The dotted curve is obtained by assuming that the fluctuations of the loop areas are small in 
each of the three peaks. An approximate value for the scaled loop area can be assigned to each of the three peaks, 
(Aq) /4:Ho fa 0.00, {Ai) /47?o ~ TOeq/2, and {A2) /4iJo ~ Woq- With these choices we assume that the magnetization 
switches when the absolute value of the external field is close to a maximum. The average loop area is calculated by 
weighting these values by the fraction of events in each peak. 
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(7.6) 



1 and {A2 



where "HF" stands for "high-frequency." This expression breaks down in the low-frequency hmit, as /2 
becomes strongly frequency dependent. Next we discuss this effect. 

For any finite time series there is a sufficiently low frequency such that the magnetization always switches during 
every half-period of the field. For very low frequencies the magnetization switches, on average, before the field reaches 
its extreme value during every half-period. Therefore, the loop-area distribution becomes unimodal, and the mean 
loop area decreases with decreasing frequency. Since the probability density of switching times, P{t), is narrow and 
unimodal for low frequencies (refer to Fig. we can choose the typical time when the magnetization switches during 
each half period, tg, to be either its mean, mode, or median. We use the latter because it is analytically more tractable 
and has an analytic asymptotic approximation. To determi ne t he median of P{t) we must solve F{ts)=l/2, where 
the cumulative probability distribution F{t) is given by Eq. (|4.6| ). To simplify the form of this equation we note that 
for low frequencies the hysteresis loops are square (refer to Fig. ^d)), which means that the growth time becomes 
small compared to tg and can be ignored. Thus, the equation for tg becomes 



F{U) = 1 - cxp 



p{t')dt' 



After rearranging and substituting x=^o/[HQsin-{ujt)]'^ ^, we obtain 



(7.7) 



ln2= / p{t')dt' 
Jo 



Po- 



JO 



[iJosin(wi')] e [«o )i -i- 1 rfi' 



= Po- 



K + d 

X ''-I e" 



zdx 



(7.8) 
(7.9) 
(7.10) 



As the frequency is dec rease d the switching field, Hs=Hq svaiiots), decreases. For very low frequencies x ^ :zq/Hq ^ 



and the radical in Eq. (7.10) is of order one, resulting in the following expression 



In 2 = po 



K+l 
'0 



Hr\d -1)0; J ^ 



K + d _ 

X e ^dx , 



(7.11) 



which may also be derived through a linear approximation to the sinusoidal field, H{t) ~ H^uit. After rearranging 
and simplifying this becomes |5G| 



CHouj = r 1 - 



K + d So 



where T{a, x) is the incomplete gamma function ^T\ , and we define 

-Soil) 



(7.12) 



c = 



In 2 (d- l)e "0 S, 
Po 



(7.13) 



With o?=2, K—Z, and the values found in Table |, C=0.101 J ^ MCSS. For small uj the hysteresis loops are practically 
square, so the scaled loop area in the low-frequency (LF) limit can be expressed as 



LF 



''cq 



Uo ■ 



(7.14) 



The switching field Hs{uj) is obtained from a numerical solution of Eq. (7.10), and the result for (^)lp /4J?o is shown 
as the solid curve in Fig. |ll|(a)-(b). Figure ^(b) shows the good agreement between this parameter-free calculation 
and the MC data for frequencies 1/R < 0.05. 
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One can obtain an approximate analytic solution by taking the first term in the asymptotic expansion |87 



T(a,x) ^ x''-^e-'''\l + - — ^ + . . .1 (7.15) 
L X i 

for large x. By ignoring the factor x°'~^, one obtains a completely analytic solution for Hs in the extreme LF limit, 
resulting in the asymptotic, logarithmic frequency dependence for the loop area 

(^)lp « 4So^ [- HCHoLu)] '~ . (7.16) 

In Fig. |ri| we show calculations of both this asymptotic analytic result (shor t-dashed curve) for the loop area and 
the loop area obtained from the numerical solution (solid curve) of Eq. (7.10). The dashed lines in Fig. ^l](b) are 



linear least-squares fits to the full numerical solution over almost four decades in frequency. The slopes of these fit 
lines would give an "effective" exponent b for the loop area, A oc u^. For the frequency regime shown, these effective 
power-law exponents for the loop area depend on frequency. These results show no evidence of an overall power-law 
relationship between the frequency and the loop area. 

To the best of our knowledge, the present study is the first to give a complete solution of A and emphasize its 
possible significance for the low-frequency power-law behavior of A reported in the literature. This issue is discussed 
in further detail in a separate paper |50(| . Theoretical considerations of nucleation effects on hysteresis have been 
mentioned previ ously in the literature ||6|,|26|,|62[ , reporting the same asymptotic frequency behavior for the loop area 



as in Eq. (7.16). However, we stress that the purely logarithmic dependence of the loop area on frequency and 



amplitude of the external field will approach the exact solution obtained from Eq. (7.12) only for extremely low 
frequencies, as illustrated by the poor agreement between the short-dashed and solid curves in Fig.|ll](b). The units 
of frequency used in our simulations may be converted roughly into Hertz by equating a phonon frequency with an 
inverse Monte Carlo step per spin, u=l MCSS~^, with i/=10^-10^^Hz as a reasonable estimate. Then, the lowest 
frequency calculated for the solid curve in Fig. 0(b) corresponds to a field period on the order of thousands to millions 
of years. The full asymptotic behavior of the loop area is realized for frequencies lower still, and is therefore well 
outside the range of feasible experiments. However, the lowest-frequency MC data point corresponds to a field period 
on the order of microseconds to seconds, suggesting the low-frequency SD behavior of the loop area described by our 
full numerical calculation could be experimentally observable. 

In extensive simulations, Acharyya and Chakrabarti have studied the hysteretic response of I sing systems with 
respect to changes in field amplitude and frequency, temperature, system size, and system dimension [p2| . In particular, 
they have presented power-law scaling relations for A. However, their simulations in d=2 are all done with very large 
field amplitudes. For the temperatures used in their papers, these field amplitudes place their system in a regime 
we refer to as the "strong-field (SF) region," where the size of a critical droplet becomes comparable to the lattice 
constant a, and the droplet picture of metastable decay breaks down [|5|. The investigation of hysteresis in kinetic 
Ising models by Lo and Pelcovits is not restricted to the SF region. In particular, the range of field amplitudes 
used places their simulations in the SD, MD and SF regions as the amplitude is increased. The hysteretic response, 
including the loop area, requires qualitatively different theoretical descriptions in each of these regions. This drastic 
change in the Ising model contrasts with a coherent rotation model of magnetization reversal, for which increasing 
the field strength would merely decrease the energy barrier that the system must overcome for the magnetization to 
be aligned with the field. For the kinetic Ising model and other spatially extended systems however, increasing the 
field not only decreases the free-energy barrier for forming a single critical droplet of the stable phase, it also changes 
the reversal mode. Our results emphasize the importance of a knowledge of the decay mode in order to obtain the 
correct frequency dependence of the average loop area over a wide range of frequencies. 

Mahato and collaborators |^,^ also considered a driven bistable system in the presence of noise. From the RTDs 
they calculate what they describe as an "average hysteresis loop," and they propose a maximum in this quantity with 
respect to noise strength as a manifestation of SR. We find their definition of a loop area rather unphysical, except for 
very low frequencies. However, the role of the hysteresis-loop area as a measure of energy dissipation indicates that a 
maximum in A may be considered an aspect of SR, which we suggest could be termed "stochastic energy resonance." 

Relatively few studies have considered B, the correlation of the magnetization with the external field. For Ising 
models subject to oscillating fields, these studies 0,^ltl have found a temperature at which the correlation attains a 
maximum value for a given frequency, amplitude, and system size. However, as recently pointed out by Acharyya [ pO[ , 
this feature is not a sign of SR. These simulations are mostly performed above Tc in a regime where the metastable 
state no longer exists. Therefore, comparison of our simulations or analytic results with these studies is not relevant. 

Our theoretical derivations of B are similar to those for A and are also given as high- and low-frequency approxi- 
mations. Figure |l^ shows our MC data for B along with the theoretical result for the high-frequency regime (dotted 
curve) and the low- frequency regime (solid curve). We assume that Hg ~ Hq for high frequencies, and that the 
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magnetization switching is stochastic with either zero, one, or two switches every period. The only non-zero con- 
tribution to the correlation comes from those periods when the magnetization switches only once. We use notation 
analogous to the high-frequency loop-area calculation for each of these three contributions, (Bq) /[2Hq/it) w 0.0, 



(Bi) /(2^/^) 
to Eq. (7.6), 



/2, and {B2) /(2Hq/t:) « 0.0, with (S)jjp /(2iJo/7r) calculated as a weighted average similar 



(^)hf 



= /o(c^)l 



+ /l(^) 




(7.17) 



Figure |lj shows that (-B)jjp / {2Hq/tt) — > as — > and as w — > 00, which is a direct consequence of the behavior of 
hiiv) (seeFig.|l](b)). 

The low-frequency calculation for B is also similar to that for A and assumes that the magnetization switches 

abruptly at |_ffs| < Hq during each half-period. This assumption together with the definition of B, Eq. ([7^), yields 



(B) 



LP 



2Ho/Tr 



Ho 



(7.18) 



As for A, the switching field Hs{uj) is obtained from a numerical solution of Eq. ( |7.10|) . There is good agreement 
with the MC data for the two lowest frequencies. As the frequency increases the theoretical result approaches zero. 
However, the assumption that the magnetization switches during every half-period begins to break down around 
l/i?=0.05, where the MC result for B passes through zero. 

Although the system studied here is both stochastic and highly nonlinear, the physical significance of the integrals 
A and B can be clarified by comparison with deterministic linear response theory. In that limit one easily finds 
that A/{ttHq) and 2B/Hq correspond to the dissipative and reactive parts of the complex linear response function, 
respectively. It is therefore natural to combine A and B into an analogous nonlinear response function. 



X{Ho,T,iu) 



1 



2B 



-A 



(7.19) 



The maximum in A and the sign change in B, which occur close together in frequency, are characteristic behaviors of 
the dissipative and a reactive parts of a response function near resonance. It is reasonable to associate this behavior 
with SR. However, as we pointed out at the end of Sec. VB, m{t) remains essentially synchronized with H(t) as the 
driving frequency is lowered further below this narrow frequency range. The system then switches reliably during 
every half-period of the field, while the switching occurs earlier and earlier in the half-period as the frequency is 
lowered. In this low-frequency regime, the norm of X remains close to its maximum value of Arric^/ {Hqtt) ^ and its 
phase gives a meaningful measure of the period-averaged phase lag. The latter increases monotonically from zero at 
w = to 7r/2 at the frequency where B crosses zero. We believe the system should be considered as resonant in this 
whole range of low frequencies. 



B. Period- Averaged Magnetization 

The period-averaged magnetization Q has been considered as a "dynamic order parameter" for systems exhibiting 
hysteresis ]3^j5^ -|57| . Those studies of the Ising model have suggested the existence of a dynamic phase transition 
between (Q) 7^ and (Q)=0. As for the hysteresis- loop areas, the statistical properties of the period-averaged 



magnetization in the SD region are not well characterized simply by its mean. Figure 13 shows the probability 
densities of Q in the SD region. For all but the very highest values of i?, the distributions show two sharp peaks near 
(5=±meq due to m{t) oscillating near the spontaneous magnetization during most of the field cycles. The contributions 
to the Q distributions near Q=Q occur when the magnetization switches twice in one period. The contributions to the 
peak centered around (5=-)-0.5 occur for those periods in which the magnetization switches only once. Note that there 
is not a corresponding peak at (5=— 0.5. This is an effect of the way we calculate the period-averaged magnetization, 
which considers the beginning of a period to start when H{t)=Q and H > 0. We would have obtained a peak near 
Q=-0.5 if we had started with H{t)=0 and H <0. 

Even by inspecting the distributions for Q, no dynamic phase transition can be seen. While the means of the 
distributions for high (low) frequencies are nonzero (zero), this happens smoothly as weight shifts from the peaks near 
Q—±l to the peak at Q—0. As we will show in a forthcoming paper, the situation is quite different in the MD region, 
where we find strong evidence for a dynamic phase transition [Bl| . 
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VIII. DISCUSSION 



The mechanism by which a metastable phase decays depends sensitively on the system size, the temperature, and the 
strength of the apphed field. For small systems and weak fields, the decay proceeds through the nucleation and growth 
of a single droplet of overturned spins. This regime has been termed the single-droplet (SD) region. In this region 
the magnetization response consists of rapid transitions between two states; one with the majority of the spins up, 
and one with the majority of the spins down. The resulting time series is well described in terms of a Poisson process 
with a time-dependent rate obtained from the nucleation rate and growth velocity of droplets of the stable phase. 
The time dependence enters the nucleation rate by replacing the constant field H by H{t)=Hosm{ujt). This central 
idea provides the analytic framework for theoretical descriptions of the quantities measured from our MC simulations. 
These quantities include residence-time distributions (RTDs), power spectral densities (PSDs), hysteresis- loop areas, 
and the correlation between the magnetization and the oscillating field. The agreement between all of our theoretical 
calculations and the MC data is very good, especially considering that the theory contains no adjustable parameters. 
All of the constants used are either known from droplet theory or are measured from MC simulations of field reversal 
in kinetic Ising models. To the best of our knowledge, the present study is the first which explicitly considers hysteresis 
for the Ising model in the SD regime. 

The frequency dependence of the RTD peak shapes and peak strengths are calculated by numerically evaluating 
analytic expressions obtained from a time-dependent extension of classical nucleation theory. The good agreement 
between our theoretical calculations and MC data supports the model of magnetization switching as a Poisson process 
with a variable rate, given by substituting a sinusoidal field dependence for the static field in the nucleation rate. 

The frequency dependence of the RTD peak strengths, the hysteresis-loop areas, and the correlation between the 
magnetization and the field all indicate the presence of stochastic resonance (SR) in the two-dimensional Ising model 
in this parameter regime. This observation is consistent with recent studies of SR in other systems of coupled bistable 
elements, some of which pointed out the importance of nucleation of kink-antikink pairs to what has been termed array 
enhanced stochastic resonance (AESR). For any nucleation process, there should be crossovers between coexistence 
(CE), single-droplet (SD), and multidroplet (MD) types of behavior as a consequence of the interplay between the 
sizes and separations of the critical fluctuation(s) and the size of the system. We believe these crossovers should 
be relevant to the dependence of the amount of enhancement on the number of elements observed in other systems 
exhibiting AESR as well. 

We also calculate the hysteresis-loop area. A, in the low- and high-frequency regimes. Because of its role as a 
measure of the energy dissipation in the system, this is a quantity of particular experimental significance. For high 
frequencies the loop-area distributions are trimodal due to the stochastic switching behavior. In this regime we 
calculate (A) as a weighted average of the loop areas obtained when the magnetization switches zero, one, or two 
time(s) during a period of the field. For the low-frequency regime we obtain an analytic expression for (A). Our 
theoretical calculation agrees well with our MC results and predicts an extremely slow crossover to a logarithmic 
dependence of the loop area on Hqlo. The switching dynamics is dominated by nucleation and indicates no overall 
power-law dependence for the loop area on field amplitude and/or frequency, in contrast to what has been claimed 
in other simulational and experimental studies. However, we emphasize that numerical analysis of data generated 
by our analytic solution, even over several frequency decades, could easily lead to the conclusion that the data were 
taken from a power law. 

The period-averaged magnetization, Q, has been proposed as an "order parameter" associated with a dynamic 
phase transition in kinetic Ising models. However, in the parameter range studied here, the probability densities of 
Q show no sign of a sharp transition as the frequency of the external field is varied. Indeed, due to the multi-peaked 
nature of the distributions for intermediate frequencies, the mean value of Q is not a useful quantity in the SD region. 
In the MD region the behavior of Q is radically different, as we will discuss in forthcoming papers 

We also computed the power spectral densities (PSDs) from the simulated magnetization time series. We qualita- 
tively explain the various features of the spectra in the full frequency range from the lowest observable frequencies 
to the rapid fluctuations due to thermal noise. Specifically, we make a connection between the RTDs and the low- 
frequency behavior in the PSDs through the characteristic time of the RTDs. Inverse characteristic times that are 
smaller (larger) than the frequency of the applied field correspond to large (small) low-frequency components in the 
PSDs. Our theoretical derivation of the characteristic time also agrees well with the characteristic times obtained 
from the simulated RTDs, except at very high frequencies. The relatively poor agreement between the MC data and 
the theory for high frequencies of the external forcing field is common to most of the quantities measured and is most 
likely due to the poor quality of the MC data for high frequencies. However, it could also be a sign of breakdown 
in the adiabatic approximation underlying the assumption that the functional form of the nucleation rate and the 
calculation of the growth-time corrections do not change for high frequencies . 

In summary, we have studied stochastic hysteresis in the kinetic Ising model, a spatially extended, bistable system 
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with thermal fluctuations. We emphasize not only the detailed differences between hysteresis in mean-field models 
and Ising models, but also the qualitatively different response that the Ising model displays for particular regimes 
of system size and field amplitude. Our theoretical and numerical study is the first to consider the effects of these 
different decay regimes on hysteresis, which may be relevant to the interpretation of simulational and experimental 
results. Especially for certain technological applications, an Ising system should be a good candidate to model the 
behavior of ferromagnetic and ferroelectric materials in oscillating external fields. Finally we note that the quantities 
that we have analyzed numerically could all be measured in experiments on hysteresis in a variety of systems and 
analyzed by methods essentially identical to our analysis of the MC data. 
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APPENDIX A: DERIVATION OF THE RESIDENCE-TIME DISTRIBUTIONS 



The a**^ residence time is defined as 



A" 



(Al) 



where the times and 9^ are shown schematically in Fig. |l^. We define i| as the time when a switching event 
takes place, as measured from the first time at which H(t)=0, after the previous switching event. Without loss of 
generality, this time can be set to t=0. 9 is the time from a switching event to the next change in the sign of the 
external field. This decomposition of the residence time facilitates the calculation of the probability density functions 
(pdf) for both and 6*". When falls during the j—1 period (see Fig ^for an explanation of the indexing scheme), 
its probability density is given by Eq. (4.14), i.e. p|(t|)=P(t|). We easily generalize to the case when i| falls during 
the jth period. This is obtained by finding the probability that, given the magnetization has not switched in the 
previous j — 1 periods, it switches at a time t^ during the jth period. 



p^{t^)^[Pno^{io)V-^P 



,27r 

t^-{j-l) — 



(A2) 



The pdf for 9, pe{9), is calculated horn p^{t^) by using the fact that 0" and 6*"+^ should be independent and identically 
distributed. Then the following substitution holds. 



t? + r+i 



ij 



1 27r 

2 CO 



(A3) 



which gives 



Pe(r)=p,(r+^) 



2 uj 



2 UJ 



oo 

p(^-r)^P„.(-r^ 



(J - 1) 



2n 



(A4) 



Thus, 



18 



(A5) 



The pdf of the total residence time, A=f| + 9, is given by the convolution of the pdfs of each term: 



n(A) 



[1 - Pnot(w)] 



(A6a) 
(A6b) 



where j = \uiA/{2tt)]. The notation \x] is defined as the smallest integer greater than x. For j=l, the integration 
limits are given by 



OaiA) = Max 
0b{A) = Max 



0,A 



0,Min[A-to,- -io] 

CO 



(A7a) 
(A7b) 



which ensure that the integrand in the equation for n(A) is positive. To implement the calculation of the RTDs, 
n(A) is numerically integrated for j=l at 50 equally spaced values in the interval < A < 27r/w, to generate the 
first peak in the RTD. To obtain the second peak, the values of the distribution are shifted to the next interval 
27r/cL) < A < Att/lj, and reduced by a factor of Pnoti^^)- This process can be repeated to obtain all of the higher-order 
peaks in the RTD. 
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Parameters 


Constants (theory) 


Constants (simulation) 1 


Ho 


O.IJ 


Ho(r) 


0.506192 J 




3.15255 


L 


64 


K 


3 (exact) 


V 


(0.465 ±0.014) J-^MCSS"' 


T 


O.STe 




(r) 


2058 MCSS 






Hosp 


(0.11 ±0.005) J 






r 


0.672 



TABLE L Parameters and constants used in this work. The values of the parameters Hq, L, and T have been selected such 
that switching occurs via the single-droplet mechanism, while the maximum nucleation rate is not too low to obtain reasonable 
simulation statistics. The constants Ho(r) and K are calculated from droplet theory [|68|-[7l[ for two-dimensional Ising systems. 
The constants (r) and r are measured from field-reversal MC simulations with the Glauber dynamic (using the parameters 
listed above). The constants fl |7^ and ly |7^] have been measured in other work (for clarity, we do not explicitly show the 
temperature dependence of these quantities in the table or elsewhere in the paper). The value for Husp is taken from Fig. 11 

of Ref. lei. 
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FIG. 1. The free energy, F{m,H,T), shown vs. magnetization m for the nearest-neighbor Ising ferromagnet on a 64 x 64 
square lattice at T=0.8Tc. Data are shown for H=Q.QJ and O.IJ. The data was obtained from a study in which a multicanonical 
MC algorithm was used to find F{m, 0, T) The inset shows an expanded view of the portion of the free energy curve near 
the metastable state for H=0.1J. The barrier height is on the order of IksT. 
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FIG. 2. Short initial segments of the magnetization time series m{t) (sohd line) and the external field H{t) (dashed line) 
vs. time t in the SD reeion. for r=0.8r.,. d—2. L—64. and Hn = O.IJ. The total length of the time series is approximately 
16.9 X ] . The time series are 

shown f 



4 
3.5 

3 
2.5 

3 

5^ 2 
i.5 
1 
0.5 




_I I I L_ 



4 





R =1.5 


ojto = 1.68 




R = 2.5 


(jjtQ = 1.23 




R =5 


oito = 0.84 




R=iO 


ojto = 0.58 




R =100 


oito = 0.18 



' I I I \-^< \ u-1- I I i"- ..\..'\-' 



"2 



3n 
4 



FIG. 3. The switching probability density, P{t)/uj vs. cut, for five values of the peri od of the external field, ii=1.5, R—2.5, 
R=5, R=10, and 7?=100. The plots are obtained from a numerical evaluation of Eq. (4.14). The inset shows the decay rate, 
p{t) vs. u!t. The time to is the earliest time during a period, for which P{t) is non-zero. Even though the decay rate always 
has a maximum when the phase equals 7r/2, the value of the phase for which P{t)/uj is maximum depends on the frequency. 
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FIG. 4. Residence-time distributions (RTDs). The time axis is scaled by the period 27r/aJ of the external field for each value 
of R, so that the peaks are centered around odd half-integer multiples. The RTDs are shown for (a) 7?= 10 (150 bins, 1239 
events), (b) R=5 (250 bins, 1439 events), and (c) R=2.5 (500 bins, 1089 events). The filled circles are obtained from the MC 
simulations. The solid curves represent the theoretical calculation presented in Appendix pU 
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FIG. 5. Peak strengths in the RTDs vs. scaled frequency 1/R. The different curves, from top to bottom, correspond to 5*1, 
<S'2, S3, a nd S4. The statistical errors are everywhere on the order of the symbol size or less. The solid curves, obtained from 
Eq. (5.4), result from the same parameter-free calculation as the RTDs. 
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FIG. 6. The inverse characteristic time of the RTDs, [27r_Rr; (t)]^^, vs. the scaled frequency, 1/ R. The inverse characteristic 
time is given in units of 10~^MCSS~^. The solid dots are calculated from MC data for the peak strengths of the RTDs. The 
inset shows examples of \nSj vs. the peak number j, along with their linear least-squares fit lines. The three frequencies shown 
in the inset are 1/R=GA, 0.2, and 0.1, denoted by diamonds, stars, and squares respectively. The error bars on the data are 
calculated from the standard deviation in the slopes of the fitted lines. The solid curve is obtained from the full numerical 
calculation of Eq. (5.7). The horizontal long-dashed line results from a low-frequency approximation obtained by setting fg=0 . 
This horizontal line is located at a value of [luRrj (t)]~^ =1.4 x 10~^MCSS~^ obtained by numerical integration of Eq. (5.11). 
The short-dashed line represents the frequency, (r) /R, which has been scaled to have the units of MCSS"'^. The open oval 
around the point for 1/_R=0.5 is a reminder of the poor statistics and large systematic error in the RTD for this frequency. 
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FIG. 7. Power spectral densities (PSDs). Spectra are shown for three different frequencies of the external field, and are 
plotted with an arbitrary offset for clarity. The inset shows the same spectra without the offset to illustrate how all three PSDs 
fall onto the thermal noise background at high frequencies. In addition to the change in the amount of smoothing, the right-hand 
section of each spectrum contains only one data point out of every 25 to facilitate plotting. The magnetization is sampled 
every 1.0 MCSS, so the Nyquist frequency is 0.5 MCSS"^ The lowest frequency that can be resolved is 2.38 x 10"^ MCSS"^ 
The dashed line with slope —2 is a guide to the eye. The arrow indicates a frequency of 1.4 x lO^^MCSS"^ in the PSD. This 
frequency value is the half-width of the spectrum predicted in Sec. V C for low frequencies of the external field. 
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FIG. 8. Probability densities for the hysteresis-loop area, A=— ^ m{H) dH. The values of R shown here are R=2, 2.5, 3, 
4, 5, 6, 10, 20, and 100. 
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FIG. 9. Representative hysteresis loops corresponding to different numbers of switching events during one period of the 
applied field. Each loop in panels (a)-(c) consists of data from a single period in a time scries with R=A. For Aq, no 
magnetization reversal occurs. This example shows a "spike" in the magnetization during one half of the period. For Ai, one 
magnetization reversal occurs during the period. For A2 in panels (c) and (d), two magnetization reversals occur. Panel (d) 
shows a hysteresis loop of type A2 from a very low-frequency time series with i?=100. Note: the arrow tails with filled dots 
indicate the start of the field cycle. 
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FIG. 10. Frequency dependence of the fraction of events contained in the peaks (a) Aq, (b) Ai, and (c) A2. The soUd 
dots are measured directly from the data in the loop-area distributions, and the dashed lines are guides to the eye. The solid 
curve in each plot is calculated using Eqs. (^a-c) with Pnot(a;) obtained from a numerical evaluation of Eq. ( |4.1^ ). (d) The 
loop-area distribution for _R=6 is shown to illustrate the peak sizes in the trimodal distribution. 



31 




32 




FIG. 11. (a). Mean and standard deviation of the loop-area distributions shown vs. the scaled frequency, 1/R. The solid 
dots are the means of the distributions shown in Fig. |^. The vertical bars are not error bars, but give the standard deviations 
of those distributions. The large standard deviations for 1/R > 0.1 indicate the trimodal distribution of the loop areas, 
typical of the SD region. For the points at the two lowest frequencies the magnetization switches nearly every half-period, 
giving a loop-area distribution which is close to unimodal. The solid and dotted curves come from two separate theoretical 
calculations. The solid curve results from numerical solution of an analytic expression for the switching field that assumes 
switching occurs during every half-period. The dotted curve comes from a calculation which uses the same values of Pnot(i^) 
used for the theoretical curve in Fig. |l^. (b). Log- log plot for the very lowest frequencies in panel (a). The solid curve is the 
same as the solid curve in panel (a). The long-dashed line segments represent linear least-squares fits to different portions of the 
numerical solution data, each covering almost four decades in frequency. The data that yield the effective exponent bi =0.077, 
are centered around log j^g(l/_R)=— 3.35; those that yield 62=0 .033 are centered around log j^q(1/-R)=— 13.78. The solid dots are 
MC simulation data. The short-dashed curve represents Eq. (7.16) with d=2 and C=0.101J~^MCSS. 
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FIG. 12. Mean and standard deviation of tlie correlation, B shown vs. the scaled frequency, 1/R. The solid dots are 
the means of the correlation distributions. The vertical bars are not error bars, but give the standard deviations of those 
distributions. The solid and dotted curves are obtained from similar calculations as those for the loop areas in Fig. hll 




FIG. 13. Probability densities for the period-averaged magnetization, Q. The values of R shown are R=2, 2.5, 3, 4, 5, 6, 
10, 20, and 100. 



34 





FIG. 14. Schematic diagrams for calculation of the RTDs. 
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